home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Languages / RLaB 1.18c / toolbox / section.r < prev    next >
Encoding:
Text File  |  1994-09-23  |  3.9 KB  |  129 lines  |  [TEXT/ttxt]
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Calculate the section properties (area, center of
  4. //              geometry, moment of inertia) of a region.
  5.  
  6. // Syntax:    section ( P , coord )
  7.  
  8. // Description:
  9.  
  10. //    The properties of polygonal sections may be calculated using this 
  11. //  rfile.  
  12. //  
  13. //  The (x,y) coordinates of the vertices of the region are stored 
  14. //  in P.  P is a two column matrix,  its first column stores xs and
  15. //  its second column stores ys.   The coordinates of the vertces must
  16. //  be stored sequentially for a complete clockwise path around the 
  17. //  region.  Holes inside the region can be subtracted by tracing 
  18. //  around them in counter clockwise direction.
  19. //  If coord is "polar", then the coordinates in P are in polar 
  20. //  coordinates (r,theta), where theta is in degree.
  21. //
  22. //      Reference: (many more)
  23. //      Wojciechowski, F., "Properties of Plane Cross Sections,"
  24. //      Machine Design, Jan. 1976, p.105.
  25. //
  26. // Example 1:  Section properties of a rectangle
  27. //
  28. //  A = [0,0;0,4;2,4;2,0;0,0];
  29. //  S = section(A);
  30. //  show_prop(A,S);
  31. //  pause();
  32. //
  33. // Example 2:  Section properties of a triangle
  34. //
  35. //  P = [0,0;0,1;1,0;0,0];
  36. //  S = section(P);
  37. //  show_prop(P,S);
  38. //  pause();
  39. //
  40. // Example 3:  Section properties of an unit circle
  41. //
  42. //  theta = (360:0:-5)';
  43. //  r = ones(theta.nr,1);
  44. //  S = section([r,theta],"polar");
  45. //  show_prop([r,theta],S,"polar");
  46. //  pause();
  47. //
  48. // Example 4:  Section properties of a C section
  49. //
  50. //  P = [1,75;ones(42,1),(80:285:5)';1.2,285;ones(42,1)*1.2,(280:75:-5)'];
  51. //  S = section(P,"polar");
  52. //  show_prop(P,S,"polar");
  53. //  pause();
  54. //
  55. // Example 5:  Section properties of a square with a center hole
  56. //  s2 = sqrt(2);
  57. //  P = [10,0;10*s2,-45;10*s2,225;10*s2,135;10*s2,45;10,0;...
  58. //       5*ones(73,1),(0:360:5)';10,0];
  59. //  S = section(P,"polar");
  60. //  show_prop(P,S,"polar");
  61. //  pause();
  62. //
  63. // Todo:
  64. //  1. make P a linked-list of several regions (may be holes), return
  65. //     the total properties.
  66. //  2. add thickness and mass density to each region in 1., so the mass, 
  67. //     the center of mass, and the mass moment of inertia can be 
  68. //     calculated.
  69. //   
  70. // Tzong-Shuoh Yang  8/10/94   (tsyang@ce.berkeley.edu)
  71. //-------------------------------------------------------------------//
  72. section = function( P, coord )
  73. {  
  74.    local(P,B,N,N1,np,prop,q);
  75.    global (eps, pi)
  76.  
  77.    if (P.nr <= 2) 
  78.    {
  79.       error("section.r: A polygon must have at least 3 points.");
  80.    }
  81.  
  82.    // if the polygon is not closed, let the last point = 1st point
  83.    if (P[P.nr;1] != P[1;1] || P[P.nr;2] != P[1;2] )
  84.    {
  85.        P = [P;P[1;1],P[1;2]];
  86.    }
  87.       
  88.    // convert polar coordinates to cartesian coordinates
  89.    if (exist(coord)) 
  90.    {  if (coord == "polar")
  91.       {
  92.          P[;2] = P[;2]*(pi/180);
  93.          P = [P[;1].*cos(P[;2]),P[;1].*sin(P[;2])];
  94.       }
  95.    }
  96.    
  97.    np = P.nr - 1;
  98.    N  = 1:np;
  99.    N1 = 2:P.nr;
  100.    prop = <<>>;
  101.    B  = P[N1;1].*P[N;2] - P[N;1].*P[N1;2];
  102.    prop.area = sum(B)/2; 
  103.    prop.cgx  = B'*(P[N;1]+P[N1;1])/6/prop.area;
  104.    prop.cgy  = B'*(P[N;2]+P[N1;2])/6/prop.area;
  105.    prop.Ixx  = B'*(P[N;2].^2 + P[N;2].*P[N1;2] + P[N1;2].^2)/12;
  106.    prop.Iyy  = B'*(P[N;1].^2 + P[N;1].*P[N1;1] + P[N1;1].^2)/12;
  107.    prop.Ixy  = B'*(P[N;1].*(2*P[N;2]+P[N1;2]) + P[N1;1].*(2*P[N1;2]+P[N;2]))/24;
  108.    // centroid values
  109.    prop.Ixx0 = prop.Ixx - prop.cgy^2*prop.area;
  110.    prop.Iyy0 = prop.Iyy - prop.cgx^2*prop.area;
  111.    prop.Ixy0 = prop.Ixy - prop.cgx*prop.cgy*prop.area;
  112.    // pricipal values
  113.    if (abs(prop.Ixy0) < eps) {
  114.       q = 0;
  115.    else if (prop.Ixx0 == prop.Iyy0) {
  116.       q = pi/4;
  117.    else
  118.       q = 0.5*atan(2*prop.Ixy0/(prop.Iyy0 - prop.Ixx0));
  119.    }}
  120.    prop.pIxx  = prop.Ixx0*cos(q)^2 + prop.Iyy0*sin(q)^2 - prop.Ixy0*sin(2*q);
  121.    prop.pIyy  = prop.Ixx0*sin(q)^2 + prop.Iyy0*cos(q)^2 + prop.Ixy0*sin(2*q);
  122.    // polar moi
  123.    prop.J0    = prop.pIxx + prop.pIyy;
  124.    // convert angle to degree
  125.    prop.angle = q*180/pi;
  126.    
  127.    return prop;
  128. };
  129.